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Abstract 

The aim of this paper is twofold. First, we set up the theory of elastic matter sources within 
the framework of general relativity in a self-contained manner. The discussion is primarily based 
on the presentation of Carter and Quintana but also includes new methods and results as well 
as some modifications that in our opinion make the theory more modern and transparent. For 
instance, the equations of motion for the matter are shown to take a neat form when expressed 
in terms of the relativistic Hadamard elasticity tensor. Using this formulation we obtain simple 
formulae for the speeds of elastic wave propagation along eigendirections of the pressure tensor. 
Secondly, we apply the theory to static spherically symmetric configurations using a specific 
equation of state and consider models either having an elastic crust or core. 


PACS: 04.40.Dg, 26.60,+c, 97.10.Cv, 97.60.Jd 


1 Introduction 

The presence of a solid crust in a neutron star is well established on theoretical grounds [I]. There are 
also speculations of solid cores due to the dominant strong interaction in the deep interiors of such stars. 
Although the presence of a solid crust is not thought to influence the mass or radius of a neutron star 
substantially there are still a number of other observational issues attributable to such a crust. The best 
observational data of what goes on inside a neutron star crust comes from pulsar timing measurements, 
occasionally showing sudden spin-ups in the otherwise very regular spin-down pattern. These events, called 
glitches, are thought to be due to interactions between the crust and superfluid neutrons penetrating the 
lattice see e.g. [2]. Another interesting feature observed in some pulsars is a smooth modulation of the timing 
measurements. This could be caused by free precession of a non-axisymmetric configuration of the crust 
with respect to the angular momentum axis. For discussions on this phenomenon see e.g. mu]. Of course, a 
crust has many other interesting implications including for instance its effects on the (quasi) normal modes 
of the star. For a nonrotating newtonian perfect fluid star all axial modes are trivial in the sense that they 
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are degenerate at zero frequency and hence correspond to non-oscillatory flows. The presence of the crust 
breaks the zero-frequency degeneracy of the axial modes as well as adding new modes due to the interface 
to the fluid core. For calculations and discussion within the Newtonian framework see 0EJ and references 
therein. Axial modes of relativistic initially isotropic elastic stars have also been considered in some detail 
in 0 ■ For some very recent relativistic calculations within the Cowling approximation and using a Hookean 
equation of state see 0. 

So far, all studies have been based on some kind of weak field approximation. Either it has been assumed 
that the gravitational field is weak or not coupled to matter oscillations, or that the deviation from an 
unsheared material state is small and can be treated linearly by a Hookean approximation. To be able to 
study nonlinear effects that fall outside such approximation schemes it is desirable to develop the general 
relativistic theory of elasticity to the point where it can be used in practical calculations. We intend for this 
work to be a step towards that goal. 

Although classical newtonian elasticity theory has a long history, going back to the 17th century and 
Hooke’s law, the development of a theory of elasticity adapted to general relativity took off rather late. The 
first demand for such a theory came in the late 50s with Weber’s bar antenna for gravitational wavesJS], 
necessitating an understanding of how these waves interact with elastic solids. For this application it is 
sufficient to use the linearised version of Einstein’s equation with a background that can be treated newto- 
nianly. Indeed, this week-field approximation was used by Weber and was further developed in subsequent 
works by DysonJTOj and others. Following unsuccessfuljTT], partially successfulpT2| and successful yet rather 
coordinate-clutteredjd attempts, a fully developed non-linear theory of elasticity adapted to general rela¬ 
tivity was given in 1973 in a paper by Carter and Quintana jT3] (hereinafter referred to as CQ), which still 
is a standard reference in the field. However, as later noted by Carter [T5|. the basic theoretical framework 
of their theory had already been given by SouriauJTB], a work that largely had passed unnoticed. Also 
preceding CQ was Maugin who has made substantial contributions to the field, including the development 
of a theory of polarizable elastic media [T71 QJ, which should be important for realistic modeling of neutron 
stars since they are believed to have strong magnetic fields. Carter also gave a discussion of such matter 
in G3- In more recent years, the general relativistic elasticity theory has been reconsidered by Magli and 
Kiiowski OH EHI who put emphasis on the gauge character of the theory, and by ChristodoulouETI. Very 
recently a number of existence and uniqueness theorems have been proved by Beig and Schmidt l22l . 

We would like to emphasise that the above bibliographical remarks by no means give a complete history 
of the subject of general relativistic elasticity theory. For further references see for instance mm 

In this paper we set up the theory of general relativistic elasticity in a self-contained manner taking as our 
starting point the theory described in CQ. Although our treatment follows CQ rather closely and sometimes 
in detail, new methods and results are presented along the way. We do not want the reader to believe that 
we consider CQ to be superior to other sound presentations of the theory in any fundamental way. We 
simply choose it as a standard reference to avoid introducing yet another completely new set of notations 
and definitions into a field which already suffers from being a bit of a jungle when it comes to comparing 
results from different workers. Also, being relativists, we favor a formulation of the theory in terms of pure 
spacetime tensor equations, rather than formulations that keep material space indices and coordinates in 
the formulae. 

The theoretical part of the paper is followed by an application to static spherically symmetric configu¬ 
rations. Having a thorough control and understanding of such models is essential as they provide the back¬ 
grounds for more elaborate applications, e.g. slowly rotating stars via the Hartle-Thorne formalism 125 ) 
or stellar oscillations of various kinds. Spherically symmetric elastic matter models appears to have been 
first studied by Magli and Kijowski^^ and were later also investigated by Park|27|. However, in this paper 
our aim is to give a recipe for numerically producing physically relevant prestressed stellar models with 
speeds of wave propagation that are explicitly known. Such models will serve as backgrounds in subsequent 
investigations of stellar perturbations. The paper is organised as follows: 

In section Q we introduce and discuss the general concepts of the theory, based on the foundation that 
the matter field is a mapping from the four-dimensional spacetime manifold to an abstract three-dimensional 
manifold. The latter will be referred to as the material space, since its points may be thought of as the 
material particles. To insure that the theory is invariant under coordinate transformations on both manifolds, 


2 


the Lagrangian of the theory, which is directly given by the energy density, will be a function of scalars that 
can be formed by contractions of the inverse spacetime metric with covariant tensors that are pulled back 
from the material space. We feel that this is the most natural way to make the coordinate invariance inherent 
to the theory, as well as to have a clear distinction between the spacetime metric and the material space 
mapping as the gravitational and matter field variables of Einstein’s equations with elastic source. This 
approach is slightly different from the one used in CQ, where the energy density is viewed as a function of 
the spatial part of the spacetime metric (which depends both on the spacetime metric itself as well as on the 
matter field), implicitly or explicitly contracted with contravariant tensors that can be identified as material 
space tensors. Using the spatial metric as a “variable” in this way would make it cumbersome to fix the 
spacetime metric and only consider variations of the matter field, which has been a source of critique from 
other workers (c/. the appendix of 22 and also | 22 ]|). However the critique is partially unjustified since it 
seems to be based on the assumption that the spatial metric is to be viewed as the elastic matter field. 

Furthermore, the basic mathematical tools that we use are standard tools of differential geometry such as 
the pull-back and push-forward associated with the material space mapping. This means that we are using 
a somewhat lighter mathematical machinery than CQ who are explicitly making use of the isomorphism of 
the spacetime tangent subspaces orthogonal to the material flowlines and the tangent spaces of the material 
space manifold, as well as an extended version of the convected derivative first introduced by Oldroyd J3J. 
The reason that we are not discussing or using these concepts is not because we by any means feel that they 
are irrelevant, but because we, during the progress of this work, found that they were not needed for any of 
our purposes. For a thorough treatment of them, we instead refer the reader to CQ. The important results 
of this section is the general form of the stress-energy tensor and its specialization to the case when a fixed 
metric is the only material space tensor field which is used to set up the equation of state, an assumption 
which will be made throughout the paper unless explicitly stated otherwise. 

In section 03 we discuss the matter equations of motion, i.e. the divergence-free condition for the stress- 
energy tensor. To this end we introduce a spatially projected connection on spacetime which can be viewed 
as the “pull-back” of the Levi-Cevita connection associated with the material space metric. This projected 
connection is given by a three index tensor field which we call the relativistic elasticity difference tensor. 
We then show that the Euler equations, to which the matter equations of motion reduces, can be put in an 
elegant form involving the four index relativistic Hadamard elasticity tensor and the relativistic elasticity 
difference tensor. To our knowledge, this form of the Euler equations is first given here. We also discuss 
how to generalise this method to more general equations of state involving material space tensor fields other 
than a metric. 

In section 0] we derive various useful relations that involve the eigenvectors and eigenvalues of the pull¬ 
back of the material space metric. This leads us to introduce three scalars that we refer to as the linear 
particle densities which simply are defined as the square roots of the eigenvalues. The advantage of working 
with these linear particle densities is that we arrive at several formulae that closely resemble analogous 
formulae for perfect fluids when the latter are described in terms of the (volume) particle density. 

In section [S] we derive the characteristic equation that determines the speed of wave propagation in 
relativistic elastic media. We do so by employing the standard Hadamard method of considering field 
discontinuities across a hypersurface, in a manner similar to Carter[29j. However, by starting out from our 
form of the nonlinear equations of motion we obtain a more direct link between the full nonlinear theory 
and the theory of wave propagation. The characteristic equation is explicitly solved in the case when the 
wave propagates in a principal (i.e. eigenvector) direction and we obtain a strikingly simple formula for the 
speed of the longitudinal waves, closely analogous to the formula for the speed of sound in perfect fluids. For 
the transversal waves, which have no analogue for perfect fluids, the propagation speed is found to be given 
by a simple algebraic formula involving the principal pressures and the linear particle densities. Although 
wave propagation in prestressed nonlinear elastic solids has been treated relativistically also by Maugin in a 
number of papers[2El 1201EH E2] j we feel that our approach and usage of the linear particle densities make 
the derivation as well as the resulting formulae much more transparent. 

In section|6]we introduce a particular class of equations of state which much resemble the quasi-hookean 
type of equations of state described in CQ, the sole difference being that we use a different definition of the 
shear scalar which we feel is more convenient. Explicit formulae for the principal pressures and speeds of 
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sound propagation are given. 

Sectional is concerned with the degenerate case when two of the linear particle densities coincide. The 
reason why we pay particular attention to this nongeneric case is that any spherically symmetric configuration 
by necessity belongs to it. 

In section[H]we set up Einstein’s equations for static spherically symmetric configurations with an elastic 
matter source. The full set of equations of motion are found to form a (nonautonomous) system of three 
first order ordinary differential equations. This should be contrasted to the perfect fluid case for which the 
number of equations is two. 

In section El we discuss the boundary conditions that static spherically symmetric configurations must 
satisfy. In particular we show that transitions between a material phase interior to an elastic solid form 
a one-parameter set of phase transitions parametrizable by the discontinuity in the tangential pressure. 
Physically speaking, the nature of such a phase transition is determined by the interaction between interior 
and exterior phases as well as the history of the configuration. However it is not directly given by the 
matching conditions of the Einstein equations. This is in contrast to the perfect fluid case where the phase 
transitions are given by the equation of state in a direct way. 

We then proceed in section m to integrate some specific toy models of neutron stars using a relativistic 
polytropic equation of state as description of the unsheared behaviour of the matter. We let the shear 
modulus be proportional to the unsheared pressure and study the effects of on the one hand including a solid 
core, and on the other having a rigid crust. As expected the equilibrium models do not differ very much 
from the corresponding perfect fluid ones for realistic values of the shear modulus. 

In the Appendix we discuss and motivate the class of equations of state introduced in sectional including 
our unorthodox shear scalar definition. The discussion is based on a second order expansion of a general 
equation of state around an arbitrarily chosen unsheared state. 

A few words on the notation and conventions used in this paper should be mentioned here. Tensors will 
be given in the abstract index notation in the manner described in Wald {TTTTi . Abstract spacetime indices will 
be denoted by lowercase Latin letters (a, 6,...), whereas material space indices will be denoted by capital 
Latin letters (A, B ,...). Greek letters (p, v ,...) or (A, E,...) will be used as indices referring to a given basis 
and for these the Einstein summation convention will not be used. In order to avoid cluttering formulae by 
explicitly writing out arguments we will use the convention that quantities depending only on the particle 
density will be denoted by a check accent. For example, the symbol p will denote unsheared pressure, to be 
explained subsequently. The metric is taken to have signature +2 and the units are such that the speed of 
light is equal to unity and the Einstein equations take the form G a b = nT a b with n being given in arbitrary 
units. 


2 Relast icity 

The theory of elastic matter sources in general relativity that we will use is based on the presentation of Carter 
and Quintana|14t Since we expect that this theory is not generally known, we will start the discussion by 
outlining the basic ideas. To begin with, we introduce an abstract three-dimensional manifold, the material 
space , A', whose points are to be thought of as the material particles when going to the continuum limit. This 
space can in general be equipped with various types of tensor fields which give the structure of the material 
in a reference state. To begin with, the only tensor field that we introduce is, using capital Latin tensor 
indices on A, a particle density form uabc = ^[abc]- This three-form is defined in such a way that when 
integrated over some volume in A', the result is simply the number of particles contained in that volume. 

A state of matter in a material filled open four-dinrensional submanifold M' C M, with M being the full 
spacetime manifold, is given by the spacetime metric g a b restricted to M' and a C 1 (at least) mapping 

: M' —> A. (1) 

This mapping if should be thought of as the matter field of the theory, analogous to a scalar field but with 
the difference that if maps to a three-dimensional space while the scalar field maps to a one-dimensional 
one. We shall always assume that the preimage if~ 1 (p) of all points p £ A' = ip(M') is a single timelike 
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curve in M', which is to be interpreted as the flowline of the particle represented by the point p. This way 
ip gives rise to an identification of the three-dimensional manifold of particle flowlines with the material 
space submanifold X' . The push-forward ip* and pull-back ip*, taking contravariant tensors from M' to X' 
respectively covariant tensors from X' to M' , are important ingredients when setting up the theory. We shall 
use the convention that the push-forward of a contravariant spacetime tensor t a ' will be denoted by t A " 
and similarly that the pull-back of a covariant material space tensor t.A... will be denoted by In other 
words we simply replace the lower case Latin spacetime indices with capital Latin material space indices and 
vice versa. From the made assumptions about ip it follows that the kernel of ip*, when viewed as a linear 
mapping of spacetime vectors v a to material space vectors v A , is comprised of the vectors that are tangent to 
the particle flowlines. It can be proved, using the commutation relations involving pull-back, push-forward, 
tensor contraction and Lie differentiation, that given a covariant spacetime tensor field f Q ..., there exists a 
fixed material space tensor field Ia... whose pullback is t a ... if and only if the following conditions hold for 
all vector fields v a that are everywhere flowline tangential; 

• All contractions of v a with any index of t a ,,, vanishes, i.e. t a ,,, is completely flowline orthogonal. 

• t a ... is Lie dragged by v a ; C v t a ... = 0. 

The pulled back particle density form n a b c = ip* tiabc is of fundamental importance for the theory and is 
used to define the flowline tangential particle current 

_ 1 -abcd^ /o\ 

^ — 3! £ Tlbcdi 

where e a i, c d is the spacetime volume form associated with g a b . By construction, n a is conserved, 

V a n a = 0, (3) 

which follows from the fact that exterior differentiation commutes with the pull-back operation, implying 
that n a b c is a closed spacetime three-form due to uabc trivially being a closed material space three-form. 
Noting that n a by assumption is timelike, we set 

n a = nu a , n = \J—n a n a , u a u a = — 1 , ( 4 ) 

where n is the particle density and u a the matter four-velocity. To proceed we introduce the spatial volume 
form 

C abc — tabcdn • (h) 

From eq. 0 it follows that 

nabc — n e a bc, ( 6 ) 

which justifies the interpretation of n as the particle density. 

An equation of state is supplied by giving the rest frame energy density p as a function of scalars that 
are formed by contracting the inverse metric g ab with pulled back covariant material space tensor fields. 
The particle density may be used as one of the scalars, as indeed, according to eq. ©, it is given by the 
contraction 



n 2 = ^ n abc n abc = g ad g be g cf n a bcn d ef ■ 

( 7 ) 

Once an 

equation of state is given, the stress-energy tensor takes the form 



Tab — PQab ^ Q gab — P^a^b “1“ Pabi 

(8) 

where 

dp 

Pab — 2 phabi ^ Pab — 0 ? 

( 9 ) 

and 

hah = UaUb H - Pab- 

( 10 ) 
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The above equations can be taken as the definition of the stress-energy tensor for elastic matter, but it is 
worth noting that T a b can also be straightforwardly derived according to 


Tab 


2 SSm 
V~ 9 s 9 ab 


( 11 ) 


from the matter action 

Sm = - J dx 4 ^/^gp. (12) 

It will prove convenient to rewrite the energy density as p = ne where e is the energy per particle. For 
instance the pressure tensor p a b is then given by the simple formula 


Pab = 2 n 


de 

dg ab ’ 


( 13 ) 


where we have used that 


dn 

dg ab 


2 Tlhab^ 


(14) 


which follows, for instance, from eq. Q. Now, scalars that are formed on spacetime by contracting g ab with 
pulled back tensors t a ,,, = ip*tA... can alternatively be formed on the material space by making the analogous 
contractions of the pushed forward tensor g AB = t^g ab = i(>*h ab with the tensors tA...- Hence, at all points 
in X we may view the energy density p as a function of g AB and fixed material space tensors tA...- Taking 
this view we shall make the assumtion that e has a minimum value e under variations of g AB such that the 
particle density n is held fixed. Such a state will be referred to as an unsheared state. This allows us to 
introduce an n-dependent tensor p AB on the material space such that g AG pcB = & A b when e = e. In other 
words, we define pab such that the minimum of e(g AB ) occurs at g AB = p~ lAB . If e, which is a function 
of n only, has a minimum at a certain particle density n = no, it follows that e has an absolute minimum 
under all variations of g AB . This minimizing state is referred to as a completely unstrained or completely 
relaxed state. We shall not assume the existence of such a state simply because the material structure in the 
interiors of a neutron star may owe its existence to the high pressures there. 

It is fairly straightforward to show that the pull-back of the Levi-Cevita volume form of pab coincides 
with the spatial volume form e a b c ■ Hence we let the volume form of pab be denoted cabc as we have 
’4 > * e ABC = e 0 be- Clearly, the relation between cabc and the particle density form uabc is analogous to eq. 
0 for their pull-backs e a bc and n a bc , i.e. 


nABC = ne A BG- (15) 

It is worth keeping in mind that, despite the notation, the particle density form uabc is in fact a fixed 
material space tensor, independent of n. It is convenient to define a new tensor conformal to pab and 
having babc as its volume form. It follows from eq. CE 3 that this tensor is given by 

k A B = n 2 / 3 p A B (16) 

Denoting differentiation with respect to n by a prime we have 

npAB = -§PAB + TAB, (17) 

where 

tab = n 1,3 k' AB (18) 

Now, since k~ 1AB k' AB = 0 due to the fact that the determinant of kAB is independent of n (since uabc is), 
we find that the trace of eq. 1171) with respect to p~ 1AB = n 2 ^ 3 k~ 1AB is 

np~ 1AB p' A B = -2 (19) 
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Thus tab , called the compressional distortion tensor , is clearly the trace-free part of nr)' AB and therefore 
measures how the deformation of the unsheared crystal structure deviates from pure conformal rescaling 
when varying the particle density. It now follows from the definition Cl that Uab is independent of n 
precisely when the material deforms conformally. In this paper we will mainly consider such materials, 
which include materials with isotropic but also, for instance, cubic unsheared structures. Consequently, we 
will have a fixed n-independent metric tensor field kAB on X. We shall also make a further restricting 
assumption, namely that the energy per particle e is a function only of the independent invariants that can 
be formed from the flowline orthogonal mixed tensor 


k a b = g ac k cb , kab = ip*kAB- 


( 20 ) 


When acting with d/dg ab on quantities that depend only on k a b 


we may reexpress this operator as 


which follows from 


d , d 

dgab ~ k <a dkb)e > 

( 21 ) 

dk c d 

9gab = d (akb)d- 

( 22 ) 


Since k a b is restricted to be orthogonal to it°, it has three independent scalar invariants, which for instance 
could be taken to be 


h = k a . 


h = k a b k b a , h = k\k b c k c a . 


(23) 


With habg being the volume form of /cab, it follows that the particle density n is in fact also a scalar 
invariant of k a b , given by 

n 2 = det ( k\) = i (I] 3 - 3 Ji/ 2 + 2J 3 ). (24) 

Rather than using {Ii, I 2 , 13 } as the set of three independent scalar invariants of k a b , it is both convenient 
and instructive to instead take the set to be comprised of n and two independent scalar invariants of the 
mixed tensor 

rf b = h ac r) cb , y ab = ip*r]AB , (25) 


whose deviation from h a b gives all information about deformations other than conformal compressions. From 
the relation k a b = n 2 ^ 3 r] a b it is clear that det (y a b ) = 1, whence ri a b has two independent invariants rather 
than three, due to its invariance under conformal compressions h ab —> C~ 2 h ab under which k a b —> C 2 k a b 
and n —> C 3 n. We shall now obtain an alternative expression for the pressure tensor by using that the 
decomposition k a b = n 2 ^ 3 ri a b of k a b into a conformal factor n 2 ^ 3 times a unit determinant mixed tensor r] a b 
gives rise to the following splitting of d/dg ab valid when acting on invariants (scalar or tensorial) of k a b : 


d i d d 

g^b = 2 nh - b Q^ + 


(26) 


where the angle brackets < , > around an index pair here denotes the symmetric and traceless part of that 
pair, where “traceless part” refers to h ab and not the full spacetime metric g a b■ Explicitly, for a tensor t ab , 
we have 

t<iab~> — i'^ab) 3 ^ c^ab- ( 27 ) 

In order to obtain eq. (El we have used eq. <d and 


dy c d 

dg ab 


h <a,y b >d 


(28) 


which follows from rj c d = n 2 ^ 3 k c d and eqs. m and El- So, by specifying e as a function of n and rj a b , 
the pressure tensor can be written as 


Pab P k ab T 7 V ab , 7T a — 0, 


(29) 
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where 


(30) 

(31) 


2 de 



This clearly displays that the dependence of e on the particle density n is directly related to the isotropic 


pressure p , while the dependence on the unimodular tensor r] a b is directly related to the trace-free anisotropy 
pressure tensor 7r ab . In particular, we see that the matter source can be interpreted as a perfect fluid iff 
the energy density is a function of the number density alone, the equation of state then being p = ne, 


p = n 2 de/dn where e coincides with e, the minimal energy per particle at fixed particle density. 


3 Equations of motion for elastic matter 


Let us now turn our attention to the equations of motion X b T ab = 0, which can also be obtained by varying 
the matter action m with respect to the material space mapping ip. For any stress-energy tensor of the 
form 



the divergence free condition splits into its projection onto and orthogonal to u a . In terms of the spatially 
projected connection D a defined by its action on an arbitrary tensor f b ' c ... according to 



the equations of motion takes the form 


P + {ph ab + p ab )Q a b — 0, 
(ph ab +p ab )u b + D b p ab = 0, 

where, as usual, a dot denotes covariant derivative along u a and 

®ab = 


(34) 

(35) 


(36) 


For elastic matter, when the stress-energy tensor is given by an equation of state according to eqs. © and 
©. it follows that eq. m is in fact automatically satisfied, leaving the u a orthogonal Euler equations as 


the only nontrivial matter equations of motions. 

When u a is hypersurface orthogonal the projection tensor h ab is the metric induced on these hypersurfaces. 
In this case the operator D a is the Levi-Cevita connection associated with h ab - In general, however, D a 
cannot be viewed as a connection defined on a family of spacetime submanifolds, but it can nevertheless be 
viewed as the Levi-Cevita “pseudo-connection” of h ab in the sense that the property 


(37) 


D a h bc — 0 


always holds. Since we are assuming that the material space is equipped with a fixed metric Uab, its 
pullback k ab provides an alternative metric on the tangent subspaces orthogonal to u a . Furthermore, we 
may also carry over the Levi-Cevita connection Da of Uab into a differential operator D a acting on spacetime 
tensors. To do so we declare that D a , like D a , should be a projected torsion-free connection with respect to 
the projection tensor h a b , i.e. there should exist a torsion-free connection V a on M such that the action of 
D a on an arbitrary tensor held f b "' c ... is given by 


D a t b -c...=h d a h b e ---h f c ---V d t e - f .... 


(38) 


It follows that D a is uniquely determined by D a and a completely howline orthogonal tensor held C c ab such 
that for all spacetime vector helds X c , we have 


D a X c = D a X c + C c ab X b . 
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(39) 



We will refer to the tensor field C c ab as the relativistic elasticity difference tensor , since the tensor field that 
gives the difference between two connections is sometimes called the difference tensor in the mathematical 
literature pH). To relate D a to D A we impose that for any spacetime vector X a and vector field Y a , it should 
hold true that 

fa{X b D b Y a ) = X b D b Y a , X B = faX b , Y a = faY a . (40) 

This fixes C c ab to be the tensor field 

C c ab = \k~ l cd (D a k bd + D b k ad - D d k ab ), (41) 

with k~ lcd being the flowline orthogonal inverse of k ab , 

k~ lac k cb = h a b , k- lab u b = k~ lba u b = 0. (42) 

Clearly, the tensor C c ab has the symmetry 


C C /TC 

ab — ^ ( ab ) • 

Moreover, as one might guess, the projected connection D a satisfies 

fa(X b D b t a -) = X B D B t A ", X B = faX b , t A - = fat a ■ 
fa(D B t A ...) = D b ta..., t a ... = fat A .... 

In particular, eq. lTl5l) implies 

Dckah — 0 . 

Letting t b " c ... be a tensor field built up out of g ab and k ab , it follows from eq. I® that 


(43) 

(44) 

(45) 

(46) 


D t b - 

c... 

implying in turn that 


dt b 


c ■■■ r> „de 


% 


de 


D a g = 


dt b - 


dg' 


de 


-{D a -D a )g de = 2 


dt b 


c ide 


_ /'Td 

dgde ° “> 


(47) 


.dt 


‘ - DJ - - < D - * ^ = 2 + • • <«) 


In particular, it follows that if the energy per particle e is a function only of scalar invariants of k a b , the 
spatial divergence of the pressure tensor that occurs in Euler’s equations m can be reexpressed as 

D bP ab = A ab cd C cd b , (49) 


where A ab cd is the relativistic Hadamard elasticity tensor 

dn ab 

A ab cd = 2-£- i -p ab h cd -h a c p b d . 

The first two terms on the right hand side of the above equation is the relativistic elasticity tensor , 


dn ab 

rpab o ^ _a6 j 

P* cd — " c\ cc [ P hcdi 


dg 


which, using eq. o. can be written in the more compact form 


E ab 


cd 


2 n 


d 

dgCd 



(50) 


(51) 


(52) 
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(53) 


With all indices raised, it can also be expressed in the familiar form (c/. CQ) 

d 2 e 


E abcd = 4 n 




which follows from eq. m and 


d ac bd, d 

= -gg 


dg ab a a dg cd ' ^ 

Both A abcd and E abcd are flowline orthogonal tensors symmetric under interchange of the index pairs ab and 
cd , i.e. 

j^abcd _ j^cdab j^abcd _ cdab 


In addition, the elasticity tensor has the symmetry 

j^jabcd _ jjj(ab)(cd) 


(55) 

(56) 


The bottom line of this discussion is that the equations of motion are given by the Euler equations (1351) 
alone, which, using eq. m . can be expressed as 

(, ph ab + p ab )ii b + A abcd C cbd = 0. (57) 


This fact will be used in section 0 when we examine wave propagation in elastic media. 

Before closing this section we would like to make a digression and argue that the Euler equations in the 
form m, which to our knowledge has not appeared before in the literature, is valid under more general 
circumstances than were assumed when deriving it above. The assumption more restrictive than necessary 
was that the energy per particle e is a function only of scalar invariants of k a b = g ac k cb with k cb being the 
pull-back of a fixed (n-independent) material space metric kcB- Indeed, the derivation of eq. (1571) can be 
immediately generalised to the case when there exists a torsion-free connection Da on X such that 

D A t B ... = 0, (58) 


for all fixed covariant tensors t A ... on X which contribute to the equation of state. This connection can be 
carried over into a projected connection D a on M in exactly the same manner as before, i.e. through eqs. 
EHt - (EH- Hence D a will again be uniquely determined in terms of D a and an elasticity difference tensor 
C c ab which, due to the requirement that D A be torsion-free, has the index symmetry (p|) . Since eqs. EU 
and E5> will hold also in this case, the pull-back of eq. E3) is of course 

D a t b ... = 0, tb... = (59) 

which implies that the remaining steps that lead up to the Euler equations in the form E3> are still valid. 

The most physically relevant situation when this generalization applies is when the reference structure 
of the material space is completely given by three crystal axis vectors 

E a a , A =1,2,3, (60) 


which, mathematically, are three linearly independent vector fields assumed to be commuting, 

[E a ,Ev} a =0. (61) 

Physically, the crystal axis vectors are to be thought of as displacement vectors pointing from one lattice 
point to three of its nearest neighbouring lattice points, with the three neighbours chosen in an ordered 
way such that the crystal axis vectors become smooth vector fields when going to the continuum limit. The 
two-dimensional surfaces that, due to the commutation requirement G3J, are guaranteed to be formed by 
any pair of crystal axis vectors (E A A , E ?, A ) with A ^ E are consequently to be thought of as crystal planes. 
The basis of one-form fields u> A A that are dual to the crystal lattice vectors, i.e. satisfying 

lo a a Ex A = 6\, (62) 
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will here be called the reciprocal lattice one-forms, as they are in direct correspondence with the reciprocal 
lattice vectors of classical lattice/crystal theory (cf. E3) where, of course, no distinction is made between 
vectors and one-forms. However, it should be stressed that a factor 27r should be inserted on the right hand 
side of eq. ED to conform to the standard definition of reciprocal lattice vectors. Now, if each lattice point is 
associated with a number Nf of particles belonging to the species that one wishes to consider as fundamental 
( e.g. baryons), the particle density form is simply 

it-abc = 3! Nf co 1 ^ of 1 b w 3 c]- (63) 

Indeed, this three-form associates a volume Nf to a cell spanned by the crystal axis vectors, 

E 1 A E 2 B E 3 c n A BC = N f , (64) 

and, clearly, such cells are in one-to-one correspondence with the lattice points. 

If we assume that no further tensor fields are introduced on X , the energy per particle e can only depend 
on the scalars 

jAS _ gab^A^B^ w A a _ - 0 * w a Aj ( 65 ) 

of which six are independent since J AS = / SA . We thus would like to introduce a torsion-free connection 

D a for which all three reciprocal lattice one-forms are covariantly constant, 

D a co A b = 0 . ( 66 ) 

This is possible because eqs. (ED and ED imply that these one-forms are closed, thus making the anti¬ 
symmetric part of eq. El identically satisfied. Consequently, eq. El implies that D A is the unique flat 
torsion-free connection for which the coordinate frame connection coefficients vanishes when using coordi¬ 
nates x A such that ui A A = d A x A . The elasticity difference tensor appearing in the Euler equations ED will 
in this case take the very simple form 


C c ab = ^ E A c D a 


w A b , 


(67) 


A—1 


where the spacetime vector fields E A a can be identified with the material space vector fields E A A since they 
are defined by 

co A aEx a =5 A v, u a E A a = 0. (68) 


if*E A a = E a a . 


(69) 


Indeed, they of course satisfy 
To arrive at eq. lilTI) we have used that 

D [a w A b] = 0, (70) 

which follows from the antisymmetric part of eq. El and guarantees that the symmetry property El holds. 
The result that the matter equations of motion in this case reduces to eq. ED with C c a b given by eq. ED 
is independent of whether or not the unsheared crystal structure is conformally rescaled under variation 
of the particle density, i.e. whether or not the compressional distortion tensor t A b is identically vanishing. 
However, if we do assume that t A b = 0, we have a preferred fixed (n-independent) material space metric 
which by necessity is of the form 


k A B = ^ As 
A,E=1 




where the metric components k A s are constants subject to the condition 

det (fc A s) = N f 2 , 


(71) 


(72) 
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thus making uabc the associated volume form of Uab- The reason why Rab can only be of this form is 
simply because there is no other way to form a fixed metric on X without introducing further tensor fields 
contributing to the equation of state. The constants /cae give the lattice structure of the material in its 
unsheared state. For instance, the special case 


fc AS = iV/^AE (73) 

corresponds to cubic lattices. The three independent scalar invariants that can be formed from k a b will 
obviously be functions of the six independent scalar invariants / As . If e does not depend on any combinations 
of the latter that cannot be formed out of k a b, the equation of state belongs to the class that this paper is 
mainly concerned with. Incidentally, this serves as a good motivation for taking the material space metric 
Uab to be flat. 

4 Eigenvalue and eigenvector formulation 

In this section we shall explore the fact that whichever scalar invariants of k a b we are using to parametrise 
the equation of state by, they can always be expressed in terms of the eigenvalues of k a b and, consequently, 
so can the energy per particle e. Since k a b is flowline orthogonal and since h a b (viewed as a three-dimensional 
metric) is positive definite, the eigenvectors of k a b corresponding to distinct eigenvalues are automatically 
orthogonal, while eigenvectors with the same eigenvalue can be chosen orthogonal. Moreover, since k a b is 
also positive definite, all eigenvalues of k a b will be positive and we shall write them as n^, p = 1,2,3, in 
terms of their positive square roots which, in a loose sense, represent one-dimensional particle densities 
in the principal directions, i.e. along the eigenvectors. We hence refer to these quantities as principal linear 
particle densities , or linear particle densities for short. According to eq. <EU), the particle density n is simply 
the product of the linear particle densities, i.e. 


n = n\n 2 U 3 . 


(74) 


Using an orthonormal basis {e^ a } of eigenvectors of k a b , implying that 


3 

Qab — U a Ui) + ^ ^ a^^ibi 
H=1 


3 

kab — ^ ^ ^ aCfibi 

fi=l 


(75) 

(76) 


it follows from eq. m that the operator d/dg ab , when acting on scalar invariants of k a b, takes the partic¬ 
ularly simple form 


9 _ i \ ' d 

Q g ab ~ 2 2^ e ^ ae ^ bn ^ dn ^ 


(77) 


The pressure tensor then becomes 


3 

Pab ^ ^ Pr- b > 

11=1 


(78) 


where 


de 


P l i=nn ll 


dn / 


(79) 


Hence p a b has the same eigenvectors as k a b, while the eigenvalues which we call the principal pressures , 
are given by eq. 
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We now proceed to calculate the elasticity tensor according to eq. m- To do so we need to know how 
to differentiate the eigenvectors with respect to g ab . Therefore we let d/dg ab act on the relations 


g c d,e p Gfj — 5pu , 

kcd&p Cf 7 Tl p 5ptj. 


This results in, respectively, 


de„ 


de„ 


pc dg ab + eac An ab e p(a e crb) 


dg a 


de„ 


2 „ ut -p 


n P e P° Q g ab +Tl °- eac PiSab ~ n P e pa.epbSpa 


dg a 


(80) 

(81) 


(82) 

(83) 


Assuming that n p ^ n a when p ^ cr, we can solve for all projections of the differentiated basis vectors. The 
result reads 


which implies 


de n 


dg a 


def 

1 „ 

Lpc dg ab 

2 e p(aCpb) 

de p c 

n 2 

,L P 

° ac dg ab 

n 2 - n£ 

— - e c e 
~ 2 

e P b + £ 


€ p (a&a b) P / 


a^p U P H 


2 ^p(a^a b) • 


Making use of eq. G3 when d/dg ab acts on the scalars p p , we find 


E ab Cd = 

P= 1 


2p f I 


d 9e P c 


P l P Q g ab 


dg 


ab 


£ 

(T—l 


dp P 


d 


Pp I ^cra^crb^p &p 


(84) 

(85) 

( 86 ) 

(87) 


Using eq. for the derivatives of the basis vectors and raising the indices a and b, we obtain our final 
form for the elasticity tensor: 


j^jabcd _ ^ ^ 

P= 1 


2 Pp e P a e p b e p c e p d + 2 £ + £ ( 

<7#p "P <7=1 V 


dp p 


na dn„ Pp) e ° e - e P e P 


Note that the symmetry property E abcd = E cdab requires that the quantities 

dpp 

- Pp 

dn a 


( 88 ) 

(89) 


(90) 


are symmetric under p <-> a. This is indeed the case, since, according to eq. G3, they can be expressed as 

dp p d 2 e de 

n a - - p p = nn p n a -— - \~nn p -— 

an„ on p on a on p 

The Hadamard elasticity tensor is now 

3(3 


dp p a b c d 

°dn a ° CT p p 
p— 1 I < 7=1 < 7 #p 


j^abcd = £ £ 


2 £ 


n P P P ~ n ° P ° r (»„ ») P (c- d)_ 


2 _ 2 e p v "e (T ''ep v “e C r“ / p p e a a e p d e p (b e a c) 


(91) 

Let us now return to the decomposition k a b = n 2 ^ 3 r] a b , which, as discussed in section|21 in a natural way 
splits Pab into its trace part ph ab and its trace-free part ir ab . Doing so, we do no longer wish to parametrise 


13 






















the equation of state by the linear particle densities but rather by the particle density n and variables 
that are closely related to the eigenvalues of ij a b- Clearly, 77 “j, has the same orthonormal eigenvectors e M “ as 
k a b, the eigenvalues being (a-j 2 , af, &£), say, with 


Tlfj, 

“ n 1 / 3 


( n l 

y?V+i?V+2 




0:10:2^3 == 1- 


(92) 


where 


n v +1 
n n +2 


cyt-i 

a H +2 ’ 


ZlZ 2 Z 3 = 1, 


and the cyclic rule p + 3 = p is being used. We find it useful to parametrise scalar invariants of rf a b 
Zfls. From n = n\n 2 n 3 and = n^+ 1 /n ^+2 it follows that 


(93) 
by the 


H / L rx ^ + 2 a 

onn on oz ^+2 


d d 

Z n+1- 


■dz, 




which via eqs. pG| and (|77| gives 

d 


Vc<c 


drf >> l 


= sE 




M =1 




M +2 


VI 


dz , 


/i+1 


It hence follows from eq. EH) that the anisotropy pressure tensor is 


7t ab — ^ fj,&fi L 1 bi 

M =1 

de 


= n 


de 


Zfi +2 


dz, 


/x+2 


' VH 


dz, 


11+1 


(94) 


(95) 


(96) 

(97) 


where the eigenvalues 717 , obey ]E /J= j 7 r M = 0 and are related to the principal pressures p M and the isotropic 
pressure p by 

Pm=P + 7V (98) 

Note that the combinations 2:^+2 d/dz fl +2 — z^+i d/dz^+i all annihilate Z 1 Z 2 Z 3 , which means that the con¬ 
dition Z 1 Z 2 Z 3 = 1 may be used to simplify anything these operators act on. 


5 Wave propagation in elastic materials 

Our discussion of the speed of sound in elastic materials will be based on the discontinuity formalism 
introduced by Hadamard about a century ago. The notation will follow Carter’s as presented in 1291 . but as 
our starting point we will take the Euler equations expressed in terms of the Hadamard tensor, see eq. ED- 
We consider a hypersurface E representing a wave front across which second order derivatives of the 
material space mapping ip are allowed to be discontinuous, whereas continuity is assumed for all first order 
derivatives of ip as well as of the spacetime metric. In other words the assumption is that ip and g a b are 
both C 1 . The metric may be assumed to have a higher degree of differentiability, but this is irrelevant to 
the discussion. The continuity of the first order derivatives of ip implies that all tensors on M which are 
pulled back from X are continuous. Since it thus follows that the particle four-velocity u a is continuous, i.e. 
well-defined on E, the wave front propagates at a well-defined speed v with respect to the particle flow, see 

fig- ID 

In order to derive an equation governing the propagation of the wave front we start by noting that the 
tensors ph ab +p ab and A abcd occurring in Euler’s equations are continuous, whereas the acceleration ii a and 
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the elasticity difference tensor C c ab may be discontinuous across E. Denoting the discontinuity of a quantity 
by surrounding square brackets we therefore have 


(ph ab +p ab )[u b ]+A abcd [C cbd }= 0. 


(99) 


Following Carter, we set 


[^a] — ^ — 1 1 


( 100 ) 


where a is the amplitude and b a is the polarization vector of the wave front. Since the discontinuity of the 
gradient of a scalar is always normal to the wave front and the Christoffel symbols are continuous due to the 
continuity of the first order metric derivatives it follows that for any tensor defined across E, there 

exists a tensor T a - b ... defined on E such that 


[V c t 5...] — b...A c , Xa — vu ai 

where v a is the propagation vector of the wave front, assumed to be normalised, 

v a v a = 1. 

It follows that there are tensors n a and ( ab such that 

[V 6 u a ] = K a X b 

[V c /c a f,] = (( a b X c • 

By contracting dim with u b and using eq. mm we find 

«“ = V 

V 

Moreover, the discontinuity of C u k a b = u c \/ c k a b + 2& c ( a Vfr)ii c = 0 gives 


a 


Cab — 2 2 ^ kc(a^b) 


which in turn implies 

[C C ab}=-^L C V a Vb- 

v z 

Finally, substituting eqs. ifTnil and JToTI) into eq. m we arrive at 

[v 2 (ph ac +p ac )-Q ac ]L c = 0, 

where Q ac = Q ( ac ) is the relativistic Fresnel tensor , given by 


r\ ac _ Aabcd. 

Lj — A v b Vd- 


( 101 ) 

( 102 ) 

(103) 

(104) 

(105) 

(106) 

(107) 

(108) 
(109) 


We have now shown that the Euler equations constrain the quantities u a , L a and v algebraically according 
to eq. « ■ In realistic situations the combination ph ab + p ab , viewed as a three-dimensional tensor, is 
expected to be non-degenerate. Clearly, the propagation vector v a can then be specified freely since for 
each v a eq. jug turns into a three-dimensional eigenvalue problem, with v 2 being the eigenvalue and i a the 
eigenvector. Since eq. (ITosl) is a three-dimensional equation, it has three independent solutions. Moreover 
it may be noted that as a consequence of the symmetry of Q ab , it follows that any two polarization vectors 
(1 )t a and ^L a corresponding to distinct eigenvalues ^v 2 and ^v 2 are orthogonal with respect to ph ab +p ab . 
For a general propagation vector v a , the Fresnel tensor can be written as 


or = E (E e ' l ° e ’ c>+ E ■ 

p=l 1 a+p 


IgiPo -Pp) 
n 2 - n 2 


L'p‘e a e a + VpV<re p K e a 


U > p dn 


= E 

p=i 
3 

+E E v p v ° 

P= 1 


®Pp , \' 

Pin . + 2^ ‘ 


<7^p 


2 n p(Pp-P<?) 

<J 9 9 

n p ~ 


dp P n 2 (p a -p p ) 

12,7 dn a n 2 - n 2 


-p ^p 


e ( a e c ) 


( 110 ) 
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where v p = e p a v a . We shall only solve the characteristic equation explicitly for the case when the propagation 
vector v a coincides with one of the eigenvectors, say e pa - The Fresnel tensor then reduces to 


Q ac =n 


dpu 


n*(p<T-Pn) 


^dn u 11 M 


er^/i 


( 111 ) 


Since 


3 

ph ac + p ac = y~^ y (p + Pa)e cr a e C r C , 

(T—l 


( 112 ) 


it follows that the characteristic equation ifTUHl) has one solution corresponding to a purely longitudinally 
polarised mode, given by 


= V u\\ = (P+Pu) 


-1„ d pr 


u »\\ ~ tP-T-W "-^ dn > 




The remaining two independent solutions are 


V 2 = = (P + Pv) 


-i n-viPv-P^) 


nl — n p 


, i a = e v a , i a v a = 0. 


(113) 


(114) 


These modes correspond to transversal waves whose polarization vector coincides with one of two eigenvectors 
e v a orthogonal to the propagation eigenvector e p a . Remarkably, the formula J33! for longitudinal waves in 
eigenvector directions exactly mimics the corresponding formula for the speed of sound c s in a perfect fluid 
with equation of state p = p(n ), p = p(n), namely 


C S 2 = (P + P ) ln (H5) 

For perfect fluids the quantity (3 = ndp/dn can be identified with the compressibility modulus, or bulk 
modulus. Thus, since a perfect fluid can be viewed as a special case of an elastic material, we take the 
analogy further by introducing the principal compressibility moduli (3 P according to 

= (116) 

OTl^ 

in terms of which the speeds of the eigenvector directed longitudinal waves (sound waves) are given by 


^11 = 




P +Pu¬ 


li is worth noting that vjj, can also be rewritten as 


dpu 

dn,, ‘ 


' dp 
dn,. 


(117) 


(118) 


analogous to the well-known perfect fluid formula c s 2 = dp/dp. 

Our resulting expressions for the nine principal speeds v^w and v p ±v are relevant to compare with those 
obtained by Maugin in 23; since in both cases it is assumed that the equation of state is set up using 
a fixed metric as the only material space structure. We could go through some lengths to establish the 
precise connection between the two sets of formulae, but refrain from doing so since it would only lead to 
the introduction and explanation of new notation rather than to something illuminating for anyone but a 
very specialised reader. In EH) the propagation speeds are given in terms of quantities that explicitly refers 
to a particular way of expressing the pressure tensor p a b by making use of the Cayley-Hamilton theorem. In 
contrast, our recipe for calculating v p i\ and from the equation of state e = e(ni, 712 , 713 ) is very simple 
and direct which should make it favourable for practical applications. 
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Although we do not treat the case of a general propagation direction, it is worth noting that the principal 
speeds and are stationary under first order variations of the propagation vector v a around the 

eigenvector e^ a . This is easily proved by taking the first order variation of the characteristic equation fToHl) 
induced by a first order variation 5v a in i/“, 

[6{v 2 )(ph ac + p ac ) - 2 A abcd + [v 2 (ph ac + p ac ) - Q ac ] Sl c = 0. (119) 

This equation is supplemented by the restrictions 

v a 5u a = 0, i a 5i a = 0, (120) 

that follows from v a and t a being unit vectors. Requiring that 5{v 2 ) = 0 gives four linear equations for the 
three components of 6i a , namely eqs. tnn and cza . In order for the system not to be overdetermined 
these equations cannot be linearly independent. To prove that we do have linear dependence in the case 
when v a is an eigenvector of k a b we contract eq. ifTHiJi with i a which gives the condition 

B bd u b Sv d = 0, (121) 

where B b d is the symmetric tensor 

B b d = A a bcd,l‘ b = B a bcdb b Pbd• (122) 

As we have seen, taking v a to be an eigenvector of k a b leads to i a being an eigenvector as well. In this case 
Bbd is diagonal in the eigenvector basis, leading to iTEHli being automatically satisfied. Thus, the number 
of independent equations for 5v a is at most three, which proves our statement that v is stationary under 
variations of propagation direction around a principal direction. 

6 A particular equation of state 

In this work will assume an equation of state where the energy per particle depends on only one invariant 
of rj a b, namely the shear scalar 

s 2 = M [(v a a) 3 - ffbtfcrfa - 24] = 5 s M’ = sOm 1 - ^) 2 ! ( 123 ) 

A»=l 

where, as before, This definition of s 2 differs from that of Carter and Quintana, for reasons 

that are discussed in the appendix. As also discussed there, the equation of state is assumed to be of the 
form 

e = e + —s 2 , (124) 

n 

where checked symbols denotes that the quantity is unsheared, i.e. depends on n only. The quantity e is 
the unsheared energy per particle that was discussed in section[2j The quantity /2, on the other hand, is the 
shear vwdulus or modulus of rigidity. Equivalently, the equation of state can be given as 

p = p+a, (125) 

with p = ne being the unsheared energy density and er = ps 2 being the shearing energy density. 

For this class of equations of states, the pressure tensor p ab takes the form 

Pab = P h ab + TTab, (126) 

where 

, - 9 de ~ n dp 

p = p+(Cl-l)a, p = n 2 —, f2 = — — (127) 

dn p dn 

lTab = \p [( V C c) 2 V<ab> - V cd Vc< a r]b>d] ■ (128) 
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The quantity p will be referred to as the unsheared pressure. In the eigenvector basis the anisotropy pressure 
tensor ir a b is given by eq. I®, the eigenvalues being 


= A (Xti+i - Xm+s) , Xu = \{z v 2 - Z D- 
It is worth noting that the .^-dependent quantities s0 and Xu satisfy 

ds; 


O &Xp r,/ 2 , 1' 

^ “ X/i ’ *d Zlx - (s ^ 


M ' Z)- 

This facilitates the calculation of the principal compressibility moduli, which we find are given by 

0y. = P + 4a - 2+ (20 - l)7r M , 

where we have introduced the shorthand notation 


dp 


/3 = n— = 0+^fi + 


(O — 1)0 + 


- dfl 


dp 


cr. 


Here, the quantity 0 is the bulk modulus 


h dp dp 

P = n ~r = (P + P)wr, 
an dp 


(129) 


(130) 


(131) 


(132) 


(133) 


and we have used nd^l/dn = 0dtt/dp. Our identification of 0 and p. as the bulk and shear moduli is justified 
in the appendix. As seen in section 0 the principal longitudinal speeds of sound are now given by 



0p 

P+Pp 


(134) 


On the other hand, the speeds of the principal transversal waves are straightforwardly calculated to be 

2 _ £ + 2ct - 5 ( 0 ^ + £7„+7T M - 7T„) 


u H±u 


P + Pu 


(135) 


Note that in the limit of zero shear, Z/Jj —> 1, the principal speeds goes over into 


Vu = 


v H±u 


V± = 7 


_ h_ 

p+p 

p 


0+ 3 A _ dp 4 fi 

dp 3 p + p 


P + P‘ 


in agreement with Carter’s results 


(136) 

(137) 


7 The degenerate case of two equal eigenvalues and orthogonal 
eigenvectors 

In this section we shall summarise some useful relations for the degenerate case when two of the eigenvalues 
of k a b are equal, i.e. when two of the linear particle densities, say n 2 and 713 , are equal. This is the 
situation one for instance encounters in the spherically symmetric case. With this application in mind we 
will label quantities corresponding to the nondegenerate and degenerate eigenvalues by r (for radial) and t 
(for tangential) respectively. We thus set 


n t = n 2 = ns 

lab — h a b r a rb — &2a^2b + ^3a^3f>i 


77/y — 77q ^ 
Ta — a-) 
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(138) 

(139) 









with t a b being the projection tensor that projects onto the degenerate eigenvector two-plane. The pulled 
back material space metric k a b can now be written as 

kab = n 2 r a r b + n t 2 t a b • (140) 

When calculating the pressure tensor we may take the view that we are dealing with an effective two- 
parameter equation of state e e ff (n r , n t ) = e(n r ,rit,rit). We are allowed to do so because the full three- 
parameter equation of state e{n\ 1 n 2 1 nz) is invariant under permutations of the arguments, which in turn 
is due to the corresponding symmetry of the scalar invariants of k a b■ Note however that, in general, we 
cannot use the effective equation of state when calculating higher order derivatives of e, since the first order 
derivatives do not obey the permutation symmetry. In particular, one needs the full three-parameter equation 
of state to calculate the elasticity tensor, which should not be surprising since it is used to determine the 
speed of wave propagation. Using the effective equation of state the pressure tensor may be written as 


Pab = Prr a r b +Pt tab, 


where 


p r = n n r 


PP ciT 

dn r ' 


p t = \nn t 


HCeff 

dn t ' 


(141) 

(142) 

(143) 


For practical purposes it is often convenient to change the two effective equation of state parameters from 
(■ n r ,nt ) to ( n,z ), with 

n = n r n 2 , z — n r /nt . (144) 

The function 2 represents the only independent quotient z p = n^+i/n^+ 2 , since clearly z\ = 1 and 
z 2 _1 = Z 3 = z. Making this change of parameters we find 


kab = n 2/3 r/ab, Vab = z l, 3 (zr a r b + 2 1 j a b)- 


(145) 


Defining a new traceless tensor 

tab = 2r a rb T j a bi 

we may calculate 7 r a b according to the very simple formula 

7^ab Q t a b 5 

where 

\ x 

q = n t = ~2 71 > = - 2 nz ~gT 

Thus, by specifying e = e e fi(n, z), the pressure tensor is given by the relation 


(146) 

(147) 

(148) 


Pab — P tlab T Q t a bi 


(149) 


where p, as before, is the isotropic pressure 


P 


= n 


tP rf: 
dn 


(150) 


As pointed out above, we cannot calculate all speeds of wave propagation using only the effective equation 
of state. However, for special cases of the directions of propagation and polarization, it is possible. Out of 
the nine speeds v^u and v in the principal directions only five are independent, namely 


1V|| = fi|| 

V t \I = ^211 = «3|| 

V r ± = V 1±2 = fl_L3 (151) 

Vt±r = ^2_L1 = ^3_L1 
Vt±t = V2±3 = ^3_L2 
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Of these five only two actually requires the use of the full three-parameter equation of state, namely those 
for which both the propagation and polarization vectors lie in the degenerate two-plane, i.e. v t \\ and Vt±t- 
The remaining three are given by 


v 


2 

r|| 


V 


2 _ 
r_L 


Vtlr 


Pr 

R - 

dp r 

P + Pr’ 

Pr - 

on r 

n t 2 (pt - 

~ Pr) 

3 q 


n 2 

1 - z 2 

np(p r 

-Pt ) 

3 q 

n r 2 - 

-np 

z~ 2 - 


(152) 

(153) 

(154) 


Having the full three-parameter equation of state of course makes it unnecessary to be cautious about 
setting 7 i 2 = n 3 before calculating various quantities, since we may simply specialise the general formulae 
to the degenerate case, i.e. set (or take the limit) 712 = 713 afterwards. Doing this for the equation of state 
treated in section El we find 


s/ = s{ = 0 

2 2 2 1 < -1 \2 
s t =s 2 =s 3 = ±{z -z) . 


= = 
M =1 


S+ • 


(155) 

(156) 


With xi = 0 identically vanishing, we denote the only independent X/j by x (without a basis index), 

X = -X2 = X3 = U Z ~ 2 - z2 ). ( 157 ) 

The anisotropy pressure scalar q is now simply given by 

q = fl X - (158) 


To conclude this section we give the principal speeds for this equation of state, 


where, as in the general case, 


v 


2 

r|| 


V 


2 

*11 




V 


2 

£_l_r 


v?±t 


/? + 4[cr- (»- \)q) 

P + Pr 

j3 + 2[cr + (Cl — 7f)q\ 

P+Pt 
fi+\((J + q) 

P + Pt 

a + H a - g) 

P+Pr 
fl + <7 

P + Pt ’ 


P — P + 3 fl + 


(n-i)n + p^- 

dp 


(159) 
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8 Application to static spherically symmetric configurations 

In this section we shall apply the elasticity theory just described to study static spherically symmetric 
configurations, with neutron stars as the obvious astrophysical objects in mind. It is remarkable that 
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although elastic spherically symmetric configurations have been studied in the past 123123 , it seems that 
not a single solution has been published, whether exact or numerical, that can be considered as a rough 
model for a neutron star. Referring to section0and using Schwarzschild coordinates, we write the spacetime 
metric as 


where 


9ab — 

’U’a'U’b -f- r a r b T t ab . 

(161) 

-e v (dt) a 


(162) 

X (dr ) a , 

n\ , 2m 

e~ 2X = 1- 

r 

(163) 

■ 2 (dn 2 ) ab 

, dfl 2 = dd 2 +siTL 2 ed(j) 2 . 

(164) 


Generically, the stress-energy tensor will only be compatible with a spherically symmetric spacetime if the 
material space metric Izab is also spherically symmetric and, in addition, the SO(3) symmetry orbits of M' 
and X' are identified by the material space mapping tjj. Hence, we set 

kAB = r A rB + tAB, (165) 


where, analogously to the above, 

r A =e x {dr) A (166) 

tAB = ? 2 (cIQ 2 )ab, dCl 2 = dd 2 + sin 2 dd(j) 2 . (167) 


The mapping V’ is now, up to irrelevant SO(3) transformations, defined through 

f = r(r ), dCl 2 = dfl 2 . 

The pulled back material tensor k a b is given by eq. tm with 

a-a dr f 

n r = e —, n t = 
dr r 


(168) 


(169) 


It is common to assume a flat material space metric, i.e. to set A = 0. However, we will keep A unspecified 
for the time being as other choices may be of interest as wellf5J|- The particle density can be written 


n = n r n 2 = ( f/ryz 

where, as in section 0 z is the dimensionless quotient 


n r _ a-a r dr 


z = — = e 
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Regardless of the matter sources, the Einstein equations for any static spherically symmetric spacetime can 
be formulated as (c/. EH) 


dv 

dr 

dm 

dr 

dp r 
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m + \nr 3 p r 
r(r — 2 m) 
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= ~(Pr + P) 


m + 1 nr 3 p r 
r(r — 2 m) 


y, d = \{Pt-Pr) 


(172) 

(173) 

(174) 


However, unless a one-parameter equation of state gives p and pt as known functions of p ri this system is 
underdetermined. In our case we are, as discussed in section 0 effectively dealing with a two-parameter 
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equation of state. An additional equation of motion can readily be found if we, to begin with, use n r and 
n t as the two independent matter variables. Indeed, eq. im directly implies that the evolution equation 
for nt is 


dn t 

dr 


1 

r 


(e x ~ x n r - n t J . 


(175) 


An evolution equation for n r can now be found by using 


which we solve for dn r /dr, yielding 


dp r dp r dn r 
dr dn r dr 


dp r dnt 
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dr /3 r \ dr dnt dr J 


(176) 


(177) 


with dp r /dr and dnt/dr given by eq. (11741 and eq. (II 7511 . respectively. Knowing p, p r and q as functions 
of n r and n t , eqs. cm cza and (fTTTI provides a closed system of three equations for three unknown 
variables m, n r and rq. The gravitational potential v is decoupled from this system and hence its governing 
equation cm can be integrated afterwards. 

From a physical point of view, it is more natural to give the equation of state in terms of n and z rather 
than in terms of n r and nt . To transform eqs. CEU and (FTTH) into evolution equations for n and z, we begin 
by noting that n 3 — n/z leads to 

dz /I dn 3 dnt 

dr y n dr nt dr 

where n/ 1 dnt/dr can be reexpressed as 


(178) 
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Moreover, by viewing p r as a function of (n, z) and using eq. cm we find 


dp r /3 r dn 3 dnt dp r 
dr n dr nt dr dz 


(180) 


which can be solved for dn/dr, resulting in 
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where we have used eqs. 1171 and <nm Inserting eq. into eq. (tT7H!i as well, the evolution equation 

for 2 is 
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1 dn 
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where dn/dr is to be inserted from eq. 11811) . Summarising, in the variables ( m,n,z ), Einstein’s equations 
take the form 
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where p , p r , q and (3 r are related by a two-parameter equation of state, given by e = e e ff(n, z) and 


P = ne, 


Pr = n 2 


de 

dn 


2 q, 


1 <9e <9p r 

q = -^nz—, (3 r = n—- + z-—. 

z oz On Oz 


( 186 ) 


We shall now specialise the system of equations Il83l) - itlSSl) to the particular equation of state introduced 
in section El To facilitate a comparison with the perfect fluid case, we use the relation dn/n = dp/p to 
replace the particle density n with the unsheared pressure p as one of the independent variables. We find 
that, with (m,p,z) as independent variables, Einstein’s equations can be put in the form 
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and where p and p are related to p through the equation of state. Note that, for vanishing shear modulus p, 
the equation for 0 is decoupled and the remaining system reduces to the standard perfect fluid form. This 
makes this explicitly closed system useful for comparing elastic models with the corresponding perfect fluid 
ones. 

It is interesting to evaluate the derivative dz/dr at zero shear z = 1. Using the fact that a = q = 0 at 
such a point we readily find 
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'm+±nr 3 p | o dp ( x _~ x 
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(197) 


Using a flat material space metric, i.e. setting A = 0, the right hand side is clearly negative for any reasonable 
equation of state, since A is then everywhere positively. Therefore, unless one chooses z > 1 at a phase 
transition, see below, z will always remain less than or equal to unity. Hence the radial pressure can never 
exceed the tangential pressure unless one chooses it to do so at the phase transition. One may suspect that 
such a choice of phase transition would tend to make the star radially unstable, but this will be the subject 
for later investigationPi?. 


9 Boundary conditions 

There are three different types of boundaries that one might consider; the center of the star r = 0, a 
transition between different materials in the star and the surface of the star r = R. At the center we require 
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that the spacetime is regular (elementary flat). This implies pressure isotropy p a b = ph a b and, given this, 
the conditions are the same as for a perfect fluid, namely finiteness of the energy density p and pressure p 
as well as vanishing of the mass function m. Since the pressure isotropy condition is fulfilled iff q = 0, eq. 
(TT5%1) shows that unless fi = 0, the function z must take its unsheared value 1, meaning that the shear s 2 
must vanish. In both cases p = p r = Pt = P and p = p hold at the center. Since the right hand sides of 
the evolution equations are all zero at r = 0 we are forced to move out the starting point of integration a 
small distance from the center when solving the evolution equations numerically. The initial conditions for 
the integration is thus given by the following expansions around r = 0: 

m = \np c r 3 + ... 

- - h (Pc + Pc) (pc+ 3p c ) -4p c p c 2 

p = p c — K Pc --;—\- r 

12 (/3 C + |Ac) 

n (Pc + Pc) {pc + 3 p c ) + 2>j3 c p c 2 
30 (Pc + |Ac) 

where a subscript c denotes evaluation at r = 0. 

The remaining types of boundary conditions are governed by the continuity of the first and second 
fundamental forms pro). which for all static spherically symmetric spacetimes leads to continuity of the 
Schwarzschild radial coordinate r, the mass function m (no surface layers), the gravitational potential v 
and its first derivative dv/dr. Via eq. (11721 1 . the continuity of dv/dr can be replaced by continuity of the 
radial pressure p r . Since we are using r as our radial variable, its continuity requirement is trivially imposed 
and will not be discussed further. 

At the surface of the star where we match the interior solution to a Schwarzschild vacuum solution of 
mass M the boundary conditions reduce to 

m\ r = R = M 

Pr\r=R = 0. 

An interesting observation is that, for fi = kp with k constant, eqf 

2| —3 — (1 — Cl s )k + \J% + 6(1 

~ ^ r=R ~ k{i + Ci s ) 

where f2 S = fl\ r —R. Remarkably, if the isotropic part of the equation of state is a relativistic polytrope, in 
which case Cl is constant (and actually equal to the adiabatic index, see section ll()[i , it follows that the value 
of z at the surface is known a priori in terms of the equation of state parameters f2 and k. 

Turning now to an interior boundary between two different phases of matter, the continuity of m and 
v can be imposed without difficulties. On the other hand the continuity of the radial pressure p r requires 
some thought since we are using the unsheared pressure p , rather than p r , as an unknown variable. Since 
the junction conditions have nothing directly to say about the shear variable z, there is in fact no unique 
way of making the matching. This means that the uniqueness result of ParkETl does not generalise to the 
case of elastic models with transitions between different phases of matter; if phase transitions are allowed for 
there is not a unique static spherically symmetric model for a given value of the central pressure p c . Since 
a real neutron star hardly can be viewed as one big piece of elastic matter, this complication is not of pure 
academic nature. Although a particular class of equations of state is considered here, it should be obvious 
that non-uniqueness is not peculiar to this class. However, we shall see below that the case when the outer 
matter phase is a perfect fluid is an important exception for which uniqueness of the matching does apply. 

Now, according to eq. ra, we are to match the following combination of the unknown variables p and 

P - | [(! - ^)(2 _1 - z) 2 + 2 (z~ 2 - z 2 )\ , (205) 
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where it is to be recalled that p and Cl are in general functions of p. For simplicity as well as for all practical 
purposes of this paper, we shall restrict ourselves to the situation when one of the phases is a perfect fluid 
(p identically vanishing) while the other is a general elastic phase. Since the equations are always integrated 
from the center and outwards, the values p- and Z- of p and z on the interior side of the boundary will be 
given while the exterior values p+ and z+ are to be determined. In fact, this implies that the cases when the 
elastic phase is interior to the perfect fluid phase and vice versa need to be treated separately. Before we do 
so, we note that for the perfect fluid phase, the value of 0 is completely irrelevant so in fact it is meaningless 
to even define z in such a region of the star. Now, starting with an elastic phase interior to a perfect fluid, 
the matching condition becomes 

P+ = V- ~ [(1 - Cl )Ur 1 - z-) 2 + 2 {zZ 2 - z*)\ , (206) 

which uniquely determines p+, while, as we just concluded, 2 + need not even be defined. If, on the other 
hand, the interior phase is the perfect fluid, the junction condition is 

p + -^± [(1 - n+)(2+ 1 - 2 + ) 2 + 2(2+ 2 - 2 2 )] = j5_. (207) 

This leads to a one parameter family of possible phase transitions characterised for instance by the value of 
z. (_. Setting 2 + = 1, it follows that p as well as all other pressures are continuous across the boundary . 1 This 
is clearly the simplest case to consider. It is interesting to observe that although z + = 1 implies continuity 
of p the converse is not necessarily true. Namely, setting p + = leads to two solutions for 2 +, 


2 1 2 

z + = 1 or 2 + = 


3 -n+ 
1 +12+ 


(208) 


where the second equality can only be valid if 12+ > 3. The second solution does not imply continuity of the 
other pressures but is probably unphysical due to the large value of 12 +. 


10 Example: relativistic polytropes 


In order to concretise the theory developed so far, we shall here consider a few examples of static spherically 
symmetric configurations. The first step is to specify the equation of state in terms of two functions of state 
e.g. p(p) and p{p). For simplicity we shall only consider models where the one-parameter relation between 
p and p is that of a relativistic polytrope, i. e. 


Pi P_ 

f-1 pi 



(209) 


where f = 0/p is the (constant) adiabatic or compressibility index and pi is the pressure which signals the 
transition between the classical polytropic p oc p T and the linear p oc p regime. More precisely, at p = pi, 
the polytropic part and the linear part contributes equally to p. This unsheared part of the equation of state 
follows from 


e = 


P lVQ 

f - 1 L 


(Vo n) 


r-i 


+ 1 


( 210 ) 


where Vo is some constant with the dimension of volume. The shear modulus p will be assumed to be 
piecewise proportional to the unsheared pressure, p = kp , where k is a dimensionless constant which is 
allowed to vanish. This assumption is supported by several calculations regarding the equation of state of 
neutron stars, see e.g. m for a review. For this class of equations of state the compressibility index f will 
be constant throughout the star. Moreover 12 will be piecewise constant and equal to T where it is non-zero. 
We chose the constant to be T = 5/3 corresponding to a non-relativistic Fermi gas. The transition pressure 

1 Note however that the total energy density p may still be discontinuous, but only if p is. 
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pi is chosen to be 1.5 x 10 37 dyn/cm 2 in order to give total radii in the order of 10 km. We shall also use 
the simplifying assumption that the material space metric is flat, i.e. we choose A = 0. 

We consider first the situation where we have a rigid crust afloat on a perfect fluid core. The transition 
from the fluid to the solid phase is taken to be approximately at the neutron drip pressure, p = 6 x 10 29 
dyn/cm 2 , and we let k take values in the range 0 < k < 0.1. Of course a zero value corresponds to a 
perfect fluid star and therefore gives a good reference model. The highest value 0.1 is probably unphysically 
large but is included in order to give insights to extreme cases. The phase transition is chosen such that all 
pressures are continuous across the matching surface. The first question we ask is how the maximum mass 
and stability properties change when we include the effect of a crust. For a few models with different k we 
have plotted the total mass M as a function of the total radius R. This allows us to read off the maximum 
mass attainable for these models. Moreover, at least in the case of a perfect fluid, an extremum of this 
curve correspond to a change in the stability for one radial normal mode[35|. Whether or not this is also the 
case when the matter is allowed to behave elastically is currently being investigated!^, but we anticipate 
that it will due to the following argument. For well-behaved equations of state the problem of calculating 
time-harmonic radial perturbations turn into a Sturm-Liouville problem with the squared frequency to 2 as 
the eigenvalues. Since the eigenvalues of such a problem can be ordered, with the lowest eigenvalue lu q 
corresponding to an eigenfunction without nodes, it is clear that instability sets in at points where ujq = 0. 
Allowing only perturbations which conserves the total mass of the star it is seen that only stationary points 
in the mass-radius diagram correspond to zero frequency modes, since it is only there the perturbed star is 
still in the equilibrium sequence, and hence does not oscillate. We point out however that if the equation 
of state involves discontinuities (as, for instance, is the case for our choice of shear modulus) one has to be 
careful since the perturbation problem will in general no longer be of the strict Sturm-Liouville type. For the 
types of discontinuities considered in this paper we do not expect that the general conclusion will change. 

In figureQwe see, as expected, that adding a solid crust does not affect the maximum mass significantly, 
even for very large k. However, for stars with large total radii, corresponding to low central pressures, the 
curves differ significantly. One should however be careful when interpreting this diagram since each model 
should be viewed as parametrised by the central pressure which is not apparent in fig. |2 For this reason it 
is useful to plot also the compactness M/R as a function of the central pressure, as done in fig. |3j Here it is 
clear that the compactness does not change significantly with increasing shear modulus, even for low central 
pressures. However, the obvious trend is that the compactness increases with increasing shear modulus. 

We turn, now, to the case when we instead have a solid matter phase in the core. We take the transition 
pressure to be 2 x 10 34 dyn/cm 2 . In this case the proportionality constant k is expected to be smaller than in 
the crust^J, if it is non-zero at all. We therefore choose to look at the interval 0 < k < 0.02, where again the 
largest value is expected to be unphysical. We again look at the mass-radius diagram (figure [U and see that 
the maximum mass actually decreases with increasing k in this case. However, looking at the compactness 
we note that the maximum value increases. This can be understood from the mass-radius diagram since the 
compactness can be parametrised by the angle between a ray from the origin to (r, m) and the r-axis. Clearly 
the maximum compactness is reached after the peak of the mass-radius curve, corresponding to an unstable 
model. One interesting possibility for neutron stars is that they, if compact enough, could allow trapped 
w-modes^ni- Based on perfect fluid calculations however, it does not seem likely that this is the case for 
real neutron st,ars[44j. Speculatively, adding a rigid core could provide the extra compactness needed. As 
discussed by Rosquist^S], it is not the total compactness that is important but rather some partial measure 
of compactness that determines if a stellar configuration allow for trapped ic-modes. A necessary (but not 
sufficient) condition is that v'r = 1 somewhere in the star. The models here does not allow for trapped 
w-modes, but a higher k in the core does indeed pushes v'r towards one. 

In order to assess the impact of the choice of parameters we have calculated a large number of models 
with various parameter values. The effects of a non-zero shear modulus do not appear to depend strongly 
on the choice of p\ or p t . Also, adding even rather large discontinuities in the tangential pressure in the 
fluid-solid phase transition does not seem to have any major impact either. 
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11 Comments and outlook 


As far as our study of static spherically symmetric solutions to the Einstein equations with elastic matter 
source is concerned, we see it as a first foundation laying step towards a study of more generic configurations. 
The obtained numerical solutions may not be highly interesting in themselves, as they, at least for realistic 
choices of the shear modulus /2, do not differ drastically from the corresponding perfect fluid model. However, 
it is important to have a good control and understanding of the background model when perturbing a 
spherically symmetric solution into less symmetric neighbours. A feature of the presented models which we 
see as oversimplified is that the equation of state is set up using a single flat material space metric throughout 
a large region of the star, stretching over several orders of magnitude in energy density and pressure. While 
the usage of a flat material space metric is naturally suggested from first principles (see the end of section 0 
one should not expect the stellar material to be ordered over such long ranges as being considered here, which 
means that the use of a flat metric as the single reference structure is not well justified. Instead, the crust 
of a neutron star is layered where each layer has a different nuclear composition and lattice structure|41l. 
Thus, a proper neutron star model should include several phases of matter, each well, at least effectively, 
described by a flat material space metric. We believe that our formalism is well suited to deal with these 
kinds of problems as well as with perturbation calculations of various kinds. 
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Appendix A: Expansion around a locally unsheared state 


To motivate our choice of invariant as well as our equation of state it is instructive to make an expansion 
of the energy around a locally unsheared state, i.e. a state where h a b = n 0 2 ^ 3 k a b for some particle density 
n = no- Using the notation °rj a b = rj a b\n=n 0 = n 0 2 ^ 3 k a b, we introduce the strain tensor e a b as 


&ab — 2 ’Hab)- 

(211) 

In the following, we shall work exclusively with tensors having mixed indices and formulate everything as 
(three-dimensional) matrix equations in index free notation. Thus needing to raise one of the indices of e a b, 
we have two natural choices to do so, either by using g ab or °rj~ lab = n 2 J 3 k~ lab . We shall for the moment 
work with both, as the former is more natural from a spacetime point of view, while the latter connects 
more closely to the classical theory and is the one used by Carter and Quintana. 

Now, with these conventions we introduce e and e to be the matrices with indices 

„a 1 (ua CL, a \ za 1/CL,—la u a \ 

e 6 - 2 \ h b- V b), e b - 2 I V b- n b), 

(212) 

or in index-free form 

e = 5(1 — °w), e=|(V 1 -l). 

Note that, formally, e and e are related by 

(213) 

c e 

6 l-2e’ C_ l + 2e 

(214) 

To proceed we expand the particle density to second order in e and e respectively, 


— ~ 1 - Tre+ i(Tre ) 2 - Tre 2 
n 0 

(215) 

~l-Tre+i(Tre ) 2 + Tre 2 , 

(216) 
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where “ ~ ” means “equal up to and including quadratic terms in e”. Note that the only formal difference 
between these two expressions is the plus sign in front of Tre 2 . So, since any quantity expanded to second 
order in e will be a linear function of Tre, Tre 2 and (Tre) 2 , it follows that such an expansion will look 
formally the same when instead expressed in e, if we use eq. (I2T51) to replace all occurrences of TV e in favour 
of n/no — 1. In the following we will therefore use only e. Now, in the classical theory of elasticity, one 
usually relates thermodynamical quantities to unit volume of the completely relaxed reference state. Hence, 
if we as our reference state take the unsheared state e = 0, which however is not necessarily completely 
relaxed, it is E = n$e rather than p = ne that should be compared to the energy per unit volume of the 
classical works. Without assuming anything about the equation of state we write the expansion as 

E~E 0 +p 0 ^-1^ +lA 0 (Tre) 2 + M0 Tr(e 2 ), (217) 

where Eq , po, Ao and po are the expansion coefficients. In general, these coefficients depend on the choice 
of no- It is apparent from eq. (I2T7I that the unsheared state can only be a completely relaxed state, i.e. a 
state that minimises E and hence also e, if 

Po = 0, A o >0, po > 0. (218) 

In this case we may identify the constants Ao and po as the Lame coefficients. The constant po is also called 
the shear modulus , while the combination /3q = Ao + ^po is called the bulk modulus. However, we warn the 
reader that this definition of the bulk modulus does not generalise to the case po ^ 0, a point which we will 
return to later in this appendix. 

We may find an interpretation for po if we note that for small e the linear term dominates over the 
quadratic and higher order terms and we have approximately the thermodynamic relation 

A E « -p 0 AV, (219) 

where A E = E — Eq, AV = — (n — ng)/n with AV being the change in the volume. Thus it follows that po 
may be identified as pressure. In situations with small pressure, for instance near the surface of a star, it 
would be sufficient to use a Hookean equation of state i.e. to let the equation of state be given by eq. TO 
with the cubic and higher order terms simply dropped and with the conditions satisfied. However, since 
we are interested in the interior of neutron stars where pressures can be extreme, a more general approach 
is called for. As noted previously in section[2]we want the scalar invariants of k used to set up the equation 
of state to consist of the particle density n along with up to two scalar invariants of the unit determinant 
matrix 17. As the invariants Tre and Tr(e 2 ) of e used in the above expansions are not exact invariants of 77, 
we need to reformulate things a little bit. So, following CQ we introduce the constant volume shear tensor 
Sab 

Sab — Vo.b)i (220) 

which in contrast to the shear tensor e a b does not depend on the particle density no of the reference state. 
As was the case with the shear tensor, there are two natural ways to raise one index, namely by using either 
h ab or r]~ lab , giving us the matrices 

s= Kl- 77 ), §=i(77” 1 -l) ( 221 ) 

Expanding an arbitrary integer power of 77 to second order in e gives 

Tr (77”) ~ 3 + 2nV, n <E Z, (222) 

where 

s 2 ~ Tr (e 2 ) - A (Tr e) 2 = Tr { [e - § (Tr e) l] 2 }. (223) 

Note that 77 has only one independent invariant up to second order in e, which therefore also will be the case 
for s. Namely, it follows that 

Tr s ~ Tr s 2 ~ s 2 ~ Tr s ~ Tr s 2 . (224) 
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It is now instructive to reexpress the energy E as a function of n/no — 1 and s 2 . Using 12151) we have, 

E — Eo + po -1^ + ^(Ao + §Ho) ^ —-+ Ho s 2 . (225) 

Due to the relations 12241) . which imply that Trs is quadratic in e to lowest order, it follows that up second 
order in e, the scalar s 2 is expressible, but actually not uniquely expressible, in terms of linear and quadratic 
invariants of s. Any function of the form s 2 ~ a Trs + (1 — a) Trs 2 -f 6(Trs) 2 will give the same energy 
contribution up to that order. For the time being we simply assume that s 2 is any such function and postpone 
a discussion of how to make an appropriate choice of it until later. To interpret eq. (12251) physically, we note 
that the first three terms, i.e. 


Eo +po (—— l') + 5 (Ao + §Ho) (—— l') (226) 

V^o ) \n 0 ) 

gives the compression part, or “perfect fluid part”, of the equation of state, while the last term hos 2 adds to 
the equation of state a resistance against shear, i.e. against matter deformations that are not pure volume 
changes. The strength of this resistance is given by the shear modulus ho- To extend this local expansion 
of the elastic matter equation of state into an equation of state valid for all admissible particle densities n, 
without neither introducing a third invariant nor inserting a nonlinear dependence on s 2 , we simply assume 
that the energy per particle e is a first order polynomial in s 2 with coefficients of arbitrary n-dependence. 
Following our previous convention that quantities depending only on n are checked, we have 


e — eo + eis 2 . 

(227) 

Since the expansion of e to second order in e is 


e ~ e 0 (n 0 ) + f-o(n 0 )(n - n 0 ) + ^€o(n 0 )(n - n 0 ) 2 + h(n 0 )s 2 

(228) 

we can immediately see that the values of the expansion coefficients at n = no are 


Eo = n 0 e 0 

(229) 

PO = «O e 'oM 

(230) 

Ao = nleo(n 0 ) - §Ho 

(231) 

Ho = n 0 ei(n 0 ) 

(232) 

Since these relations hold for any no, we may introduce n-dependent elasticity parameters according to 

p = ne 0 


P = n 2 e o 

A = n 3 e t " — | h 

(233) 

H = ne\. 


The parameters p, p and p, will be called the unsheared energy density , the unsheared pressure and the shear 
modulus, respectively. The parameter A, whose value at p = 0 is one of the classical Lame coefficients, does 

not have a useful physical interpretation and will not be mentioned by us elsewhere, 
given by the second derivative of eo will be described by the bulk modulus 

Instead, the information 

$ = = n 3 eQ + 2p = \ + 2p + | fi. 

dn ^ 

(234) 

Let us now discuss the choice of shear invariant s 2 which so far only is defined 
The primary requirements that we which to impose are: 

up to second order in e. 
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• s 2 should locally, around h a b = °rj a b for any admissible value of no, coincide with Tr | [e — i(Tre)l]“| 
up to second order in e. 

• s 2 should be positive definite in the sense that it should be be strictly positive except at zero shear 
where it should vanish. 

• s 2 should be an exact invariant of rj. 

According to (12241) a natural choice is to simply replace the locally defined shear matrix e in (12221 by either 
of the globally defined constant volume shear matrices s or s, i.e. by either 

Tr { [s - i (TV s) 1] 2 } = i Tr { [,? - i (Tr1] 2 } = | n" 4 / 3 Tr { [k - § (Tr k) l] 2 } (235) 

or 

Tr{[s- i(Trs)l] 2 } = iTrj^ 1 - ^Tr^ 1 )!] 2 } = ± n 4 / 3 Tr { [k " 1 - ±(Trk- 4 )l] 2 } (236) 

Indeed, the latter is the choice made by Carter and Quintana. However, we will choose a different definition. 
Both of the above expressions have rather cumbersome fractional power dependences on the linear particle 
densities n M . We found this to lead to unnecessarily inconvenient formulae when using one of those in 
practice. The fractional powers are clearly due to the fact that the invariants are of second degree in rj or 
r] 1 whereas we are in a three-dimensional space. This suggests that a third degree invariant should be more 
suited from this point of view. The following choice turns out to be a suitable candidate, 


s 2 = ^ [(Tr??) 3 - Tr(7j 3 ) - 24] = ± n~ 2 [(Trk) 3 - Tr(k 3 ) - 24detk] (237) 

Remarkably, this invariant looks precisely the same expressed in jy _1 , i.e. it is symmetric under rj <-> r]^ 1 . 
It is not at all obvious that it fulfills the first two of the three requirements listed above. That the first 
requirement holds true can be proved in a straightforward manner. In order to show that our s 2 fulfills the 
second requirement of positive definiteness we simply write down the invariant in terms of the linear particle 
densities. It is then given by 


s 


2 


1 


36 (m n 2 n 3 ) 2 

1 

12 


[K 2 


+ n 3 ) 


2 \3 


1 

12 


ni 

y+i 

( Tl'l V 

n 2/ 

1 1 

\ni) 



\2 

Til 

n 2 

^ + 

n 2 

Til 
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n 2 

n 3 


n 2 


(rij 6 + rc 2 6 + n 3 6 ) - 24 (m n 2 n 3 ) 2 ] 

~T + (~T -6 

n 2 J \nij \n 3 ) 

( m _ ni Y 

\ni n 3 ) 


(238) 


Clearly, the last expression shows that our s 2 is positive except at n\ = n 2 = n 3 , i.e. at vanishing shear 
rj = 1. Let us finally consider the degenerate case of two equal particle densities. In section 0 we saw that 
setting n 2 = n 3 gives 




(239) 


where z = ni/n 2 . Had we instead chosen the definition of s 2 to be given by eq. 123511 . the corresponding 
expression would be 

s 2 = \ z 2/3 {z~ 1 - z) 2 , (240) 

or, using eq. caa, 


= u- 2 / V‘- *)’ 


(241) 


For this degenerate case it is clear that the choice (I23T1) has the simplest dependence on the particle densities. 
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Figure 1: This schematic picture shows the wave front E moving in spatial direction 
v a with speed v with respect to the flow u a at a point p. It is convenient to introduce 
the vector \ a = v a — vu a which is normal to E and gives the wave propagation speed 
as its contraction with u a . Also sketched by a dashed line is the tangent space of E 
at p. 
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Mass-radius diagram 



Figure 2: Mass-radius diagram for the models with a solid crust. 


Compactness 



Figure 3: Compactness as a function of central pressure for models with a solid 
crust. Inserted is a magnification of a portion of the graph. There it can be seen 
that the compactness increases with increasing shear modulus, although the effect is 
very small. 
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